rng(328120)

load('../../01_data/05_temp/mixed_logit_rc6_totVA_DISAD_Ndis') 



% select the best Xdraws and clean them up
    [Xdraws_best] = func_likely_draws_current(IJitidx, IJchoiceset,IJcurrentFirst, Xdraws,Xmat,thetahat,y_logistic);
    
    % redo Xdraws_best to be as long as everything else
    [Ci, IAi, ICi]=unique(IJitidx(IJchoiceset==1));
    Xdraws_unique=Xdraws_best(IAi,:);
    Xdraws_best2=Xdraws_unique(IJitidx,:);
    
% draw logit errors

error_draws=evrnd(0,1,length(IJitidx),1); 

IJipref=[Xmat Xdraws_best2]*thetahat + error_draws;

% rank positions
pref_temp = IJipref;

rank_vec = zeros(size(pref_temp));

INC = accumarray(IJitidx,ones(size(pref_temp)),[],@sum);

maxC = max(INC);

for cc=1:maxC
    maxcurrent = accumarray(IJitidx,pref_temp,[],@max);
    maxcurrent = max(maxcurrent,-1000000+1);
    rank_vec(pref_temp==maxcurrent(IJitidx)) = cc;
    pref_temp(pref_temp==maxcurrent(IJitidx)) = -1000000;
end

rank_vec = rank_vec ./ INC(IJitidx);

rank_vec = 100*(1-rank_vec);

commute_q =  quantile(Xmat(:,2),20);

commute_v = zeros(size(Xmat(:,2)));
for vv=1:20 
   commute_v(commute_v==0 & Xmat(:,2)<=commute_q(vv)) = vv; 
end
commute_v(commute_v==0 & Xmat(:,2)>commute_q(20)) = vv+1;

mean_rank_commute = accumarray(commute_v,rank_vec,[],@mean);
num_rank_commute = accumarray(commute_v,ones(size(rank_vec)),[],@sum);
commute_ceil_short = accumarray(commute_v,Xmat(:,2),[],@mean);
commute_ceil_short(num_rank_commute==0) = [];
mean_rank_commute(num_rank_commute==0) = [];

p10_rank_commute = zeros(21,1);
p90_rank_commute = zeros(21,1);
for vv=1:21
    p10_rank_commute(vv) = prctile(rank_vec(commute_v==vv),10);
    p90_rank_commute(vv) = prctile(rank_vec(commute_v==vv),90);
end
p10_rank_commute(num_rank_commute==0) = [];
p90_rank_commute(num_rank_commute==0) = [];

clf

figure(1)
scatter(commute_ceil_short(2:end),p10_rank_commute(2:end),'x','black');
hold on
scatter(commute_ceil_short(2:end),mean_rank_commute(2:end),'o','black');
hold on
scatter(commute_ceil_short(2:end),p90_rank_commute(2:end),'+','black');
xlabel('commute time')
ylabel('preference percentile')
ylim([0 100])
saveas(gcf,'../../03_output/figures/rank_commutetime.png')
saveas(gcf,'../../03_output/figures/rank_commutetime_BW.eps')

% output

output_q =  quantile(Xmat(:,4),20);

output_v = zeros(size(Xmat(:,4)));
for vv=1:20 
   output_v(output_v==0 & Xmat(:,4)<=output_q(vv)) = vv; 
end
output_v(output_v==0 & Xmat(:,4)>output_q(20)) = vv+1;

mean_rank_output = accumarray(output_v,rank_vec,[],@mean);
num_rank_output = accumarray(output_v,ones(size(rank_vec)),[],@sum);
output_ceil_short = accumarray(output_v,Xmat(:,4),[],@mean);

output_ceil_short(num_rank_output==0) = [];
mean_rank_output(num_rank_output==0) = [];

p10_rank_output = zeros(21,1);
p90_rank_output = zeros(21,1);
for vv=1:21
    p10_rank_output(vv) = prctile(rank_vec(output_v==vv),10);
    p90_rank_output(vv) = prctile(rank_vec(output_v==vv),90);
end
p10_rank_output(num_rank_output==0) = [];
p90_rank_output(num_rank_output==0) = [];

clf

figure(1)
scatter(output_ceil_short(2:end),p10_rank_output(2:end),'x','black');
hold on
scatter(output_ceil_short(2:end),mean_rank_output(2:end),'o','black');
hold on
scatter(output_ceil_short(2:end),p90_rank_output(2:end),'+','black');
xlabel('total value added')
ylabel('preference percentile')
ylim([0 100])
saveas(gcf,'../../03_output/figures/rank_output.png')
saveas(gcf,'../../03_output/figures/rank_output_BW.eps')



% N disadv

frac_disadv_q =  quantile(IJNDisadv,20);

frac_disadv_v = zeros(size(IJNDisadv));
for vv=1:20 
   frac_disadv_v(frac_disadv_v==0 & IJNDisadv<=frac_disadv_q(vv)) = vv; 
end
frac_disadv_v(frac_disadv_v==0 & IJNDisadv>frac_disadv_q(20)) = vv+1;

mean_rank_frac_disadv = accumarray(frac_disadv_v,rank_vec,[],@mean);
num_rank_frac_disadv = accumarray(frac_disadv_v,ones(size(rank_vec)),[],@sum);
frac_disadv_ceil_short = accumarray(frac_disadv_v,IJNDisadv,[],@mean);

frac_disadv_ceil_short(num_rank_frac_disadv==0) = [];
mean_rank_frac_disadv(num_rank_frac_disadv==0) = [];

p10_rank_frac_disadv = zeros(21,1);
p90_rank_frac_disadv = zeros(21,1);
for vv=1:21
    p10_rank_frac_disadv(vv) = prctile(rank_vec(frac_disadv_v==vv),10);
    p90_rank_frac_disadv(vv) = prctile(rank_vec(frac_disadv_v==vv),90);
end
p10_rank_frac_disadv(num_rank_frac_disadv==0) = [];
p90_rank_frac_disadv(num_rank_frac_disadv==0) = [];

clf

figure(1)
scatter(frac_disadv_ceil_short(2:end),p10_rank_frac_disadv(2:end),'x','black');
hold on
scatter(frac_disadv_ceil_short(2:end),mean_rank_frac_disadv(2:end),'o','black');
hold on
scatter(frac_disadv_ceil_short(2:end),p90_rank_frac_disadv(2:end),'+','black');
xlabel('number economically disadvantaged')
ylabel('preference percentile')
ylim([0 100])
saveas(gcf,'../../03_output/figures/rank_Ndisadv.png')
saveas(gcf,'../../03_output/figures/rank_Ndisadv_BW.eps')

clear

rng(328120)

load '../../01_data/05_temp/mixed_logit_principal_totVA_DISAD'    

IJX = IJX(Xmat(:,2)~=0,:);
IJjtidx = IJjtidx(Xmat(:,2)~=0,:);
Xmat = Xmat(Xmat(:,2)~=0,:);


error_draws=evrnd(0,1,size(IJX,1),1); 

IJjpref=Xmat*thetahat(1:end-1) + error_draws;

% rank positions
pref_temp = IJjpref;

rank_vec = zeros(size(pref_temp));

% IJitidx         = rawdata(:,1);
% IJjtidx_teach   = rawdata(:,2);
% IJjtidx_princ   = rawdata(:,3);
% IJapp_year_cf   = rawdata(:,4);
% IJVAmadev_cf    = rawdata(:,5);
% IJVAmean_cf     = rawdata(:,7);

INC = accumarray(IJjtidx,ones(size(pref_temp)),[],@sum);

maxC = max(INC);

for cc=1:maxC
    maxcurrent = accumarray(IJjtidx,pref_temp,[],@max);
    maxcurrent = max(maxcurrent,-1000000+1);
    rank_vec(pref_temp==maxcurrent(IJjtidx)) = cc;
    pref_temp(pref_temp==maxcurrent(IJjtidx)) = -1000000;
end

rank_vec = rank_vec ./ INC(IJjtidx);

rank_vec = 100*(1-rank_vec);

% output

outputvec = Xmat(Xmat(:,2)~=0,2);
rankvecoutput = rank_vec(Xmat(:,2)~=0);

output_q =  quantile(outputvec,20);

output_v = zeros(size(outputvec));
for vv=1:20 
   output_v(output_v==0 & outputvec<=output_q(vv)) = vv; 
end
output_v(output_v==0 & outputvec>output_q(20)) = vv+1;

mean_rank_output = accumarray(output_v,rankvecoutput,[],@mean);
num_rank_output = accumarray(output_v,ones(size(rankvecoutput)),[],@sum);
output_ceil_short = accumarray(output_v,outputvec,[],@mean);

output_ceil_short(num_rank_output==0) = [];
mean_rank_output(num_rank_output==0) = [];

p10_rank_output = zeros(21,1);
p90_rank_output = zeros(21,1);
for vv=1:21
    p10_rank_output(vv) = prctile(rank_vec(output_v==vv),10);
    p90_rank_output(vv) = prctile(rank_vec(output_v==vv),90);
end
p10_rank_output(num_rank_output==0) = [];
p90_rank_output(num_rank_output==0) = [];

clf

figure(1)
scatter(output_ceil_short(2:end),p10_rank_output(2:end),'x','black');
hold on
scatter(output_ceil_short(2:end),mean_rank_output(2:end),'o','black');
hold on
scatter(output_ceil_short(2:end),p90_rank_output(2:end),'+','black');
xlabel('total value added')
ylabel('preference percentile')
ylim([0 100])
saveas(gcf,'../../03_output/figures/principal_rank_output.png')
saveas(gcf,'../../03_output/figures/principal_rank_output_BW.eps')

